;+
; NAME:
;       RUN_MORPH_MATCHING
;
; PURPOSE:
;       Run the morphological mathcing using the Hyperion model
;
; CATEGORY:
;       EMAF
;
; CALLING SEQUENCE:
;       RUN_MORPH_MATCHING
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE  [Creates an IDL save file.]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       POSTAGE_BLOCK -- Block containing the processed GLIMPSE
;                        postage-stamp images, so they may be kept in
;                        memory instead of being read in from disk
;                        each time.
;       S_BLOCK       -- Stores the Bolocat structure so it
;                        doesn't have to be read in from disk
;                        each time the routine is called.
;       VEL_BLOCK     -- Contains the BGPS velocity structure (large)
;                        so it doesn't have to be read in over
;                        and over.
;
; MODIFICATION HISTORY:
;
;       Created:  10/18/12, TPEB -- Initial Version.
;       Modified: 11/02/12, TPEB -- Forced SAVE & RESTORE verbosity to
;                                   follow /SILENT keyword.  Added
;                                   COMMON block for Bolocat to reduce
;                                   disk reads.
;       Modified: 11/06/12, TPEB -- Made compliant with use of alpha
;                                   testing routines.
;       Modified: 11/13/12, TPEB -- Added galactocentric positions to
;                                   the IDL save file.
;       Modified: 11/20/12, TPEB -- Changed "well-constrained"
;                                   parameter (CONSTR_ML) to be where
;                                   Pchoo GE 0.78.  Above this value
;                                   the larger dml error bar is <=
;                                   1.5kpc.
;       Modified: 11/27/12, TPEB -- Changed PVEC save file location to
;                                   ./emaf_paper/save_files/
;       Modified: 12/03/12, TPEB -- Added check for GRS information --
;                                   if not yet done, then do.
;
;-

FUNCTION DPDF_ML_EBAR, level, TOTPROB=totprob
  
  ind = WHERE(totprob GE level[0], nind)
  area = (nind EQ 0) ? 0.d : TOTAL(totprob[min(ind):max(ind)])
    
  ;; print,'Returning: ',abs(area/TOTAL(totprob) - 0.6826d)
  
  RETURN, abs(area/TOTAL(totprob) - 0.6827d)
END

PRO RUN_MORPH_MATCHING, DRIVEMORPH=drivemorph, _EXTRA=_extra, $
                        SILENT=silent, UPSILON=upsilon, ALPHA=alpha, $
                        ERRVLSR=errvlsr

  COMMON POSTAGE_BLOCK, common_img, common_hdr
  COMMON S_BLOCK, s
  
  message,'Starting...',/inf
  
  ;; Parse
  drivemorph = KEYWORD_SET(drivemorph)
  verbose = ~KEYWORD_SET(silent)
  suff = (n_elements(upsilon) EQ 0) ? '' : string(upsilon,format="('_',F0.2)")
  suff += (n_elements(alpha) EQ 0)  ? '' : string(alpha,format="('_a',F0.2)")
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  d = dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART
  
  ;; Load EMAF structure
  restore,'./emaf_paper/bgps_emaf.sav',VERBOSE=verbose
  n = n_elements(emaf)
  
  ;; Check for GRS information... if no, then do.
  IF fix(total(emaf.grs)) EQ 0 THEN add_grs_info
  
  ;; Define the arrays to hold distance measurements
  dml = fltarr(n,3)
  dbar = fltarr(n,2)
  tdpdf = fltarr(n,!MW.NBINS)
  pchoo = fltarr(n)
  
  ;; Name of the pvec file
  IF n_elements(pfn) EQ 0 THEN pfn = './emaf_paper/save_files/irdc_morph_pvec.sav'
  
  ;; Check if we need to drive the morophological matching
  IF drivemorph OR ~ FILE_TEST(pfn) THEN BEGIN
     
     ;; IF COMMON BLOCK IS EMPTY, fill it...
     ;; IF n_elements(common_img) EQ 0 THEN load_common_postage
     
     driver_morph,CNUM_LIST=emaf.cnum,/SAVE,TOPHAT=0,_EXTRA=_extra,$
                  SILENT=silent,UPSILON=upsilon,ALPHA=alpha, FNSAVE=pfn
     
  ENDIF
  
  ;;Load remaining SAVE files
  restore,pfn,VERBOSE=verbose
  COMMON VEL_BLOCK, velocity
  IF n_elements(velocity) EQ 0 THEN $
     restore,'./local/bgps_velocity_struct.sav',VERBOSE=verbose
  IF n_elements(s) EQ 0 THEN s = read_bgps_csv()
  
  tandist = !MW.R0*cos(emaf.l * !dtor)/cos(emaf.b * !dtor)/1.d3
  
  message,'Starting DPDF processing...',/inf
  ;; Go through each source and do the KDIST
  FOR i=0L, n-1 DO BEGIN
     IF (i+1) MOD 100 EQ 0 THEN $
        message,'Processing DPDFs for object '+string(i+1,format="(I3)")+$
                ' of '+string(n,format="(I0)"),/inf
     
     this  = WHERE(velocity.cnum EQ emaf[i].cnum, nind)
     kprob = kdist_spectrum(emaf[i].l, emaf[i].b, velocity[this].v_std, $
                            velocity[this].densegas, CLEM=clem, ERRVLSR=errvlsr)
     
     dpdf = kprob * pvec[this].invx * prob_lonlat(s[this])
     
     ;; Analyze!
     dpdf /= max(dpdf,mlind)
     
     dml[i,0] = d[mlind]/1.d3
     
     fi = WHERE(d/1.d3 GT tandist[i], nfi)
     ni = WHERE(d/1.d3 LE tandist[i], nni)
     
     ;; P_choo
     pchoo[i] = dml[i,0] GT tandist[i] ? $
                total(dpdf[fi]) / total(dpdf) : $
                total(dpdf[ni]) / total(dpdf)
     
     ;; 68.3% error bars on dml
     parinfo = {Value:0.5d,limited:[1b,1b],limits:[0.d,1.d],$
                tnside:0,parname:'Level',step:0.05}
     functargs = {TOTPROB:dpdf}
     level = TNMIN('DPDF_ML_EBAR',FUNCTARGS=functargs,/AUTODERIVATIVE,$
                   PARINFO=parinfo,STATUS=status,ERRMSG=errmsg,/QUIET)
     ei = WHERE(dpdf GE level, nei)
     IF nei EQ 0 THEN BEGIN
        print,string(i+1,format="(I3)"),level,nei
        print,m4_stat(dpdf)
        print,m4_stat(kprob)
        print,m4_stat(pvec[this].invx)
        print,m4_stat(prob_lonlat(s[this]))
     ENDIF ELSE BEGIN
        dml[i,1] = (d[mlind] - min(d[ei])) / 1.d3
        dml[i,2] = (max(d[ei]) - d[mlind]) / 1.d3
     ENDELSE
     
     ;; Moment distances
     dpdf /= total(dpdf)
     dbar[i,0] = total(dpdf * d) / 1.d3
     dbar[i,1] = sqrt( total(dpdf * d * d)/1.d6 - dbar[i,0]*dbar[i,0] )
     tdpdf[i,*] = dpdf
     undefine,dpdf,kprob
  ENDFOR
  
  ;; Find Galactocentric Positions
  omni_lbd2rz,emaf.l,emaf.b,dml[*,0]*1.d3,rml,zml
  omni_lbd2rz,emaf.l,emaf.b,dbar[*,0]*1.d3,rbar,zbar
  
  ;; Full-width ML error bar
  fw68 = dml[*,1] + dml[*,2]
  
  ;; constr_ml = (pchoo GE 0.78) OR (abs(dml[*,0] - tandist) LE 0.5)
  constr_ml = (fw68 LE 2.2) OR (abs(dml[*,0]-tandist) LE 1.0)
  
  ;; For objects at \dtan, use \dbar instead of \dml
  ;; Define d_use as the distance to use.  Set un-constrained to -100
  duse = dml
  ind = where(~constr_ml)
  duse[ind,*] = 0.
  duse[ind,0] = -100.
  
  tan = constr_ml AND abs(dml[*,0]-tandist) LE 1.0
  far = constr_ml AND ~tan AND dml[*,0] GT tandist
  near = constr_ml AND ~tan AND ~far
  
  ind = where(tan)
  duse[ind,0] = dbar[ind,0]
  duse[ind,1] = dbar[ind,1]
  duse[ind,2] = dbar[ind,1]
  
  ;; Print out statistics for the dopey user to read
  print,'           Total Sources: ',fix(n_elements(constr_ml))
  print,'Well-Constrained Sources: ',fix(total(constr_ml))
  print,'            Near Sources: ',fix(total(near))
  print,'             Far Sources: ',fix(total(far))
  print,'   Tangent-Point Sources: ',fix(total(tan))
  
  
  
   ;; Find Galactocentric Positions
  omni_lbd2rz,emaf.l,emaf.b,dml[*,0]*1.d3,rml,zml
  omni_lbd2rz,emaf.l,emaf.b,dbar[*,0]*1.d3,rbar,zbar
  omni_lbd2rz,emaf.l,emaf.b,duse[*,0]*1.d3,ruse,zuse
 
  message,suff,/inf
  save,dml,dbar,pchoo,tdpdf,rml,zml,rbar,zbar,constr_ml,fw68,$
       duse,ruse,zuse,near,far,tan,VERBOSE=verbose,$
       FILENAME='./emaf_paper/save_files/morph_distances'+suff+'.sav'
  
  undefine,dml,dbar,pvec,pchoo,tdpdf,constr_ml,rml,zml,rbar
  undefine,zbar,ruse,zuse,duse
END
